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Abstract 

Among the methods for solving ODE-IVPs, the class of General 
Linear Methods (GLMs) is able to encompass most of them, ranging 
from Linear Multistep Formulae (LMF) to RK formulae. Moreover, it 
is possible to obtain methods able to overcome typical drawbacks of 
the previous classes of methods. For example, order barriers for stable 
LMF and the problem of order reduction for RK methods. Neverthe- 
less, these goals are usually achieved at the price of a higher compu- 
tational cost. Consequently, many efforts have been made in order to 
derive GLMs with particular features, to be exploited for their efficient 
implementation . 

In recent years, the derivation of GLMs from particular Bound- 
ary Value Methods (BVMs), namely the family of Generalized BDF 
(GBDF), has been proposed for the numerical solution of stiff ODE- 
IVPs [11]. In particular, in [8 , this approach has been recently devel- 
oped, resulting in a new family of L-stable GLMs of arbitrarily high 
order, whose theory is here completed and fully worked-out. Moreover, 
for each one of such methods, it is possible to define a corresponding 
Blended GLM which is equivalent to it from the point of view of the 
stability and order properties. These blended methods, in turn, allow 
the definition of efficient nonlinear splittings for solving the generated 
discrete problems. 

A few numerical tests, confirming the excellent potential of such 
blended methods, are also reported. 
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1 Introduction 

General Linear Methods, in the form introduced by Burrage and Butcher 
|13j , are numerical methods for solving ODE-IVPs which, in order to advance 
the integration one step, require some information from the previous step 
{external stages), along with some internal stages to be computed at the 
current step. Such methods are able to describe, as extreme cases, both RK 
methods and LMF (the latter, when used as Initial Value Methods (IVMs)) 
|14| . A GLM with r external stages and s internal stages, when applied with 
stepsize h for solving the IVP 

y' = f{t,y), y(to) = yoGM-, (1) 

is usually described as 

yN = /i(^®/„)F(yN) + (t/^/^)y[n-i], (2) 

yM = /i(B^/„)F(yW) + (y®/„)y["-y, n = l,2,..., 
where: 

• the matrices A G R^^^ U G M''^'^, B G M''^^ V G M''^^ characterize 
the method; 

• yW,i^(yW) G M*™' are the vector of the internal stages and the cor- 
responding values of the function /, respectively; 

• y["^-'^l,y["l G M^™ are the vectors of the external stages. 

Many papers have been devoted, across the years, to the construction and 
the analysis of GLMs (see, e.g., the review article [16]; see also [T5 | [T71 [THl 
da EOl El [221 ES]). Clearly, the efficient implementation of GLMs depends 
on the properties of the matrix A. In particular, we shall here consider 
implicit GLMs for which the external and the internal stages coincide, i.e., 

A = B, U = V e W'', and y^ = y["l G M^™. 
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Moreover, matrix A is nonsingular and the order of accuracy of each entry of 
y["] is equal to p (to be specified later), which is, therefore, the global order 
of the method. We mention that sometimes GLMs with approximations 
having the same order have been also called "peer methods" (see, e.g., [29j). 
The approach that we shall consider is the one introduced in [HI Chapter 11, 
Section 6] and is based on Boundary Value Methods (BVMs). In more detail, 
the family of BVMs which we shall consider is that of Generalized Backward 
Differentiation Formulae (GBDF) \10\ lllj. For the efficient solution of the 
generated discrete problems, we shall consider the blended implementation 
of such methods. This implementation, introduced in [2] for block implicit 
methods (see also [31 HI El El US])) naturally induces an efficient splitting 
procedure for the solution of the generated discrete problems. The linear 
analysis of convergence for this splitting will be done according to |7j. 

With these premises, the structure of the paper is the following: in Sec- 
tion [2] Boundary Value Methods (BVMs) are briefly sketched, along with 
their block form; in Section [3| the family of Generalized Backward Differ- 
entiation Formulae (GBDF), in the BVMs class, is considered, in order to 
obtain L-stable GLMs of arbitrarily high order; in Section H] the correspond- 
ing Blended GLMs are described; Section [5] is devoted to the implementation 
details of the methods; Section [6| contains some numerical tests; finally, Sec- 
tion [7] contains a few concluding remarks. 



2 Boundary Value Methods (BVMs) and their block 
form 

Boundary Value Methods (BVMs) are a relatively new class of numerical 
methods for ODEs based on an unconventional use of LMF. Even though 
recent developments of such methods have been also obtained (see, e.g., 
\26\ I27j). nevertheless, the main reference for such methods remains the 
book |llj . The basic idea, on which BVMs rely, can be easily explained 
by just applying a fe-step LMF, which can be described, by using the usual 
notation, as 

k k 

i=Q i=Q 

to the standard test equation 

y' = Xy, Re{X) < 0. (4) 
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If we assume we know the first G {1, . . . , A;} values of the discrete solution, 
as well as the last k — v, say 

yo,---,yu-i, and yN-k+u+i, - ■ ■ ,yN, (5) 

then, under suitable conditions for the method, if 

\ziiq)\ < \z2iq)\ < ■ ■ ■< \zk{q)\ 

are the zeros of the characteristic polynomial of the method, 

k k 

p{z) — qa{z), where q = hX, and p{z) ='^^aiz\ a{z) = ''^^ I3iz\ (6) 

i=0 1=0 

then the discrete solution is approximately given by 

Vn « c [z^{q)T , for < n < A^, (7) 

with the constant c independent of N [llj. When v = k we have the 
usual way in which LMF are used; i.e., we approximate the continuous IVP 
by means of a discrete IVP. We speak, in such a case, of an Initial Value 
Method (IVM). On the other hand, when u < k we are approximating the 
continuous IVP by means of a discrete BVP. The latter defines a Boundary 
Value Method (BVM) used with [v, k—v) -boundary conditions. This freedom 
in the choice of the number ly of the initial conditions in ([5]), allows us to 
overcome the usual Dahlquist's barriers for stable LMF (see [llj, for further 
details) . 

Nevertheless, the IVP ([T]) only provides the initial condition yo, whereas 
the remaining set of fc — 1 additional conditions in ([5]) are not known. How- 
ever, they can be retrieved implicitly, by introducing a suitable set of — 1 
additional initial methods, 

^aS^\ = /.^/3?V., J = 1,..., 1^-1, (8) 

1=0 i=0 

and k — u additional final methods, 

k k 

C^f-^yN-^ = hY, f^EfN-^, j = U + l,...,k, (9) 

i=0 i=0 

which are independent of the main method (l3|). Moreover, in such a way, 
BVMs can be also implemented in block form (block BVMs [11]), which 
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is computationally more appealing: as a matter of fact, block BVMs have 
been successfully implemented in the computational codes GAM [21] and 
GAMD [32] . In more detail, the block method, with blocksize r, is used in 
a sequential fashion by solving, at the nth step, a discrete problem in the 
form 

A <S) ImY -hB® Imf = 0, 



where 



A 



a 



(1) 



a 



(1) 



(1) 



a, 



(u^l) (u-l) 



a; 



QO 



Ql 



a 



(u-l) 



a. 



(^+1) 



a 



a 



V 4^ <' c^r I 

matrix B is similarly defined by replacing the a-s with the /3-s, 

y = ( y{n-l)r7 ?/(n-l)r+l) • • • > Unr , f = ( /(n-l)r) /(n-l)r+l) • • • > fnr , 



(fe) 



(fc) 



(10) 



and y(n-i)r) f{n-i)r ^rc known from the previous step (when n = 1, yo is the 
initial condition for problem ([I])). Block BVMs allow an easier implemen- 
tation of variable stepsize, since only the stepsize h of the current "block" 
can be varied. Several families of block BVMs have then been defined (see 
dH, for full details). 



3 Generalized Backward Differentiation Formulae 
(GBDF) and GLMs 

Let us consider the particular class of fc-step LMF having the polynomial 
(t(z) in its simplest form, i.e., 

a{z) = z^, 

where j G {0, 1, . . . , k}, with the coefficients of the corresponding polyno- 
mial p{z) (see dni)) uniquely defined by imposing that the method has the 
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maximum possible order p = k. When j = k, one obtains the well-known 
family of BDF which, however, provides 0-stable methods only up to /c = 6 
and ^-stable methods up to /c = 2. Nevertheless, when 

j = iy=\{k + 2)/2-], 

and the formula is used as a BVM with (v, k — i/)-boundary conditions, 
i.e., by fixing the values dS)) of the discrete solution, then it turns out that 
the resulting methods are stable for all values of k, and have been called 
Generalized BDF (GBDF) fTOl EH H]. By the way, when k = 1,2, one 
obtains the usual first two BDF. In view of the implementation of GBDF as 
GLMs, one may assume they know the initial values of the discrete solution. 
On the other hand, the final k — v values can be retrieved implicitly by using 
a suitable set of additional final methods Q, consisting of LMF having the 
following characteristic polynomials: 



i=0 



Si). 



j = u + l,...,k, 



with the coefficients {ap^} uniquely defined by imposing the order p = k 
conditions. If we assume, for the sake of simplicity, that a constant stepsize 
h is used, then by introducing the vectors 



/ Vn+l \ 



\ Vn+r / 



n+1 



and the r x r matrices Ai and A2 such that 
/ ao 



[^1 I A2 



a. 



ao 



V 



a, 



(fe) 



Yold 



Oik 



(v. 



'n-r+l 



\ Vn 



a 



Ok 



a 



(k) 



J 



ill) 



it turns out that the discrete problem, at a given point tn = nh, can be cast 
in matrix form as 



(^2 <^ Im)yr. 



(Al (g) Im)yold- 



(12) 
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This corresponds to a GLM of the form ([2]) with coinciding internal and 
external stages, and A = A2^, U = —A2^Ai (indeed, matrix A2 turns out 
to be nonsingular) . Such methods were called Block BVMs with Memory 
(B2VM2S) in [llj. In the GLM notation, the abscissae defining this GLM 
are Cj = i, z = 1, 2, . . . , r. However, in order to guarantee the L-stability of 
the corresponding method, it is sometimes required to change them into the 
following set of abscissae: 

Ci — i, i — 1,...,-^ 1, 

i 

= + j = 0,...,r-e, (13) 

m=0 

with positive {Cm} s.t. Cr = i- 
A first possible choice, as also suggested in [IT], is that of setting 

Cm = C^\ m = 0,...,r-i, (14) 
where C is the positive root of the polynomial 

r-£ 

p{z) = z""^^ - 1- 

m=0 

However, in general such a value of C is an irrational number, which im- 
plies that the coefficients of the corresponding GLM are irrational numbers 
as well. In order to have methods whose coefficients are always rational 
numbers, a slightly different choice for the abscissae can be made, i.e.: 

im = 2T-i+i _ I ' m = 0,...,r-i. (15) 

Clearly, when i = r we have a uniform mesh, whereas for i < r the last 
r — £ + l stepsizes are geometrically decreasing, in the case of choice (I14p . or 
approximately halved, in the case of choice ()15p . In such a case, the points 
corresponding to the abscissae {cq, . . . , C£_i, c^} are equally spaced and will 
be the ones needed for the subsequent integration step. This implies that r—£ 
zero columns must be appropriately inserted in matrix ^1 in (|lip . i.e., those 
referring to the approximations at the abscissae {q, . . . , c^-i} (according to 
pT] . such points are called auxiliary points). Consequently, such a GLM 
based on GBDF is uniquely determined by the triple of integers {k,r,i). 
Clearly, the order of the corresponding method is p = k. In addition, for 
all practical values of k, by appropriately choosing the values of {r,i), it is 
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possible to guarantee the L-stability of such methods. However, it turns out 
that, for each value of k, the couples (r, £) are not unique. Consequently, an 
additional criterion will be considered, in the next section, which is aimed 
to speed-up the iterative solution of the corresponding discrete problems, 
via the blended implementation of the methods. 

4 Blended General Linear Methods 

In the previous section, we devised a procedure for obtaining a whole class of 
implicit GLMs based on GBDF, containing L-stable methods of arbitrarily 
high order. We now consider the problem of efficiently solving the generated 
discrete problems which, at a given time step, assume the form (see ([2])), 

y-h{A^I^){ = ri, (16) 

where, according to ([TT]) and (fT2]) . 

y = ynew, i = inew, A = V = Im)y old- 

A discrete problem in the form (|16p can be efficiently solved, under suitable 
hypotheses, via the blended implementation of the method. We recall that 
the blended implementation of block implicit methods has been previously 
considered in the framework of one-step methods O El HI El El El El] and 
has been implemented in the computational codes BiM [H [31] and BiMD 
121 EI]. In order to describe the blended implementation of the GLM it 
is convenient to consider its application for solving the usual test equation 
dH). In such a case, in fact, the discrete problem reduces to a linear system 
of dimension r: 

{I-qA)y = rj, q = hX, (17) 

where, hereafter, / denotes the identity matrix of dimension r. Such a linear 
system is clearly equivalent to 

^(A-^ - ql)y = jA-^n = rj^, (18) 

with 7 > a free parameter. By introducing the weight function 

e{q) = {l-^q)-'l, (19) 

we can then obtain an equivalent discrete problem by combining, with 
weights 9{q) and / — 0{q), respectively, the linear systems (fT7|) -(fT8 ]) : 

Miq)y ^ {e{q){I-qA)+ jiI-e{q)){A-'-qI))y 

= 0{q)rj + {I-6{q))ri, ^ r^{q). (20) 
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Equation (|20|) defines the Blended General Linear Method (Blended GLM) 
corresponding to the original GLM (|17p . Then, in view of the fact that 



M{q) 



/, for g ~ 0, 
-^ql, for \q\ » 1, 



the fohowing blended iteration for solving the discrete problem (I20p is nat- 
urally induced: 

iV(g)y(*+i) = {I- 7g/)y(^+^) = {N{q) - M{q))y^'^ + ri{q), ^ = 0, 1, . . . . 

(21) 

Remark 1 By considering that in [21]} N{q)^^ = 9{q) (see ( fl^) ). it is quite 
straightforward to realize that, in the case of problem the corresponding 
blended iteration formally becomes 

<5« = N-^ {e ((/ - 7^-1) ® Im (y('^ - r?) - h{A - 7/) /^f») 

+7 (y(^) -n)-hl(^l^ f ) , (22) 
y(m) ^ yW_^W, i = 0,l,..., 
w/iere 

iV = 0-i = l0(/„-/i7J), (23) 

with J the Jacobian of f evaluated at the last known point. Consequently, 
only the factorization of one matrix of the size of the continuous problem is 
required. 

Coming back to the iteration (12 ip . we observe that (see [Si |V|) the it- 
eration matrix corresponding to the blended iteration (f2T]) turns out to be 
given by 

Z{q) = / - N{q)-^M{q) = ^^^^—^A~\A - ^if. (24) 

According to the linear analysis described in we shall now study the 
spectral properties of Z{q)^ in order to obtain a linear analysis of convergence 
for the iteration (1^ . For this purpose, hereafter let p{q) denote the spectral 
radius of Z{q). Clearly, the iteration will be convergent if and only if p[q) < 
1; moreover, the set 

r = G C : p{q) < 1} 

is the region of convergence of the iteration. For the sake of completeness 
(see [7] for details), we recall that the iteration is: 
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• 0-convergent, if p(0) = 0; 

• A-convergent, if C T; 

• L-convergent with index Uoo, if it is A-convergent, 

Zoo = lim Z{q) = O, 

and i/qo is the index of nilpotency of Z^o- 

Clearly, the iteration (|2ip is 0-convergent, since ^(0) = O. Moreover, one 
easily verifies (from (I24p ) that 

— O, as g — )• oo, 

so that ^-convergence and L-convergence (with index 1) are equivalent. In 
addition to this, 




^ pq, for g ~ 0, 

w PodQ'^, for \q\ > 1, 



where p and poo are the nonstiJJ amplification factor and the stiff convergence 
factor of the iteration, respectively. Finally, for a 0-convergent iteration, 
A-convergence is equivalent to requiring that the maximum amplification 
factor, 

p* = max p{ix), 

(with i denoting, as usual, the imaginary unit) is not greater than 1. Accord- 
ing to [7\ , p, Poo , and p* are evaluation parameters for the blended iteration 
(|2ip (the smaller they are, the better the properties of the iteration). Conse- 
quently, if p* < 1, the iteration is L-convergent and, therefore, appropriate 
for L-stable methods [7j, like the GLMs defined in Section |3l 

Because of the form ()24p of the iteration matrix, it turns out that (see 
also [7]): 

p = p {A'HA - ll?) , Poo = 7"'P, P* = (27)-'p, 

with p (A~^(A — 7/)-^) denoting the spectral radius of that matrix. Con- 
sequently, the value of the free positive parameter 7 is chosen in order to 
minimize p*. In Table [H the relevant figures are reported for selected GLMs 
based on GBDF, each characterized by the corresponding triple {k,r,i), 
when the choice (|14p for the auxiliary points is considered. Similarly, in 
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Table [21 the same results, obtained with the choice ()15p of the auxiliary 
points, are listed. In both tables, we list r — i (i.e., the number of the auxil- 
iary points) in place of i. The specified values of the blocksize r, have been 
here chosen as small as possible, in order to minimize the computational 
cost per step. 

As it can be seen, the blended iteration corresponding to each method 
turns out to be L-convergent. Moreover, the value of the parameter 7 (with 
the only exception of that corresponding to /c = 4, for both the choices of 
the auxiliary steps) turns out to coincide with the value 



7 = mm 

suggested in [3] for a class of Blended Implicit Methods (see also [7]). 

As an example, we list the matrices A and U which define the third 
order GLM corresponding to the triple (3, 2, 2), requiring no auxiliary steps 
(hereafter, let c S M'' denote the vector of the abscissae). 



(1,2: 



A 



1 

23 



22 
36 



U 



1 

23 



28 
27 



(25) 



and those defining the fourth order GLM corresponding to the triple (4, 4, 3), 
with the choice (II 4|) . 



A = 

( 0.85795933248329 
1.14013555838905 
1.16454145194449 
V 1.16304178987375 



-0.22588594615620 
0.75295315385401 
0.91413597782316 
0.90057211551936 



0.09292560797716 
-0.30975202659054 
0.27338586108581 
0.50490826434742 



-0.02384741384935 
0.07949137949784 

-0.07015876367984 
0.09507592794253 



( 0.07149661104027 
0.09501129653242 
0.09704512099537 
V 0.09692014915615 



-0.44184164162566 

-0.52719452791448 

-0.53021970356702 

-0.53024220062923 



1.37034503058538 
1.43218323138206 
1.43317458257165 
1.43332205147308 



and with the choice {[T5 
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c = (l.2.*.3) , 



A 



6336684 



1 



/ 5429268 -1381941 570807 -174960 \ 

7249176 4606470 -1902690 583200 

7421568 5690784 2388672 -732160 

V 7415388 5637357 3628233 198936 / 



U = 



6336684 



1 



/ 452439 -2798388 8682633 \ 

604098 -3345408 9077994 

618464 -3365888 9084108 

V 617949 -3366036 9084771 / 



which sHghtly differ from each other. 

Finally, in Figures [T] and [21 we plot the boundary of the stability regions 
of the methods listed in Table [U whereas in Figures [3] and [U we plot the 
boundary of the stability regions of the methods listed in Table [2j 

5 Implementation details 

In the actual implementation of the above Blended GLMs, three main points 
need to be clarified: 

• the choice of a suitable starting procedure; 

• an efficient local error estimate; 

• the variable stepsize implementation of the methods. 

All of them are briefly sketched here. 

Concerning the first issue, namely the definition of a starting procedure 
to obtain the first vector of approximations, from the initial condition yo of 
problem ([1]), a natural candidate is the block GBDF of order k and blocksize 
r = k. In more detail, at the very beginning, one solves the discrete problem 



hB Imf = 0, 



(26) 



where 



y = ( yo, yi 



• • • , Vkf 



f = ( /o, /l, ■ ■ ■ , fk 



12 



k 


r 


r-l 


7 


P 


Poo 


P* 


3 


2 





0.7223 


0.2272 


0.4355 


0.1573 


4 


4 


1 


0.6195 


0.3802 


0.9908 


0.3069 


6 


5 


1 


0.6063 


0.5734 


1.5600 


0.4729 


8 


6 


1 


0.5769 


0.6380 


1.9170 


0.5530 


10 


7 


1 


0.5502 


0.6626 


2.1887 


0.6021 


12 


9 


2 


0.5271 


0.7345 


2.6438 


0.6968 


14 


10 


2 


0.5127 


0.7366 


2.8022 


0.7183 


16 


11 


2 


0.4999 


0.7345 


2.9393 


0.7347 



Table 1: Parameters of the blended iteration associated with Blended GLMs, 
based on GBDF, characterized by the triple (/c, r,tj and the choice ([1] 



k 


r 


r-(. 


7 


P 


Poo 


P* 


3 


2 





0.7223 


0.2272 


0.4355 


0.1573 


4 


4 


1 


0.6249 


0.3827 


0.9801 


0.3062 


6 


5 


1 


0.6082 


0.5740 


1.5520 


0.4719 


8 


6 


1 


0.5778 


0.6381 


1.9113 


0.5522 


10 


7 


1 


0.5507 


0.6625 


2.1845 


0.6015 


12 


9 


2 


0.5274 


0.7345 


2.6407 


0.6964 


14 


10 


2 


0.5130 


0.7366 


2.7998 


0.7180 


16 


11 


2 


0.5000 


0.7345 


2.9374 


0.7344 



Table 2: Parameters of the blended iteration associated with Blended GLMs, 
based on GBDF, characterized by the triple {k,r,£) and the choice (fT5|) . 



and, by denoting with {a^^} the as coefficients of the additional initial and 
final methods ([8]) and ([9]), the matrices A and B are defined as (see (fTOll ): 



A 



a 



a 



\ (fc) 



a 



(1) 



a 



ai 



a 



a 



a 



(1) 



a 



(^-+1) 



a 



(k) 



B 



:o|4) e 



For the efficient solution of equation (l26j) . a corresponding blended methods, 
with the associated blended iteration, can be conveniently and easily defined. 
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The estimate of the local error is obtained by deferred correction, as 
described in [TTl Chapter 10] (see also [HE]), by first "plugging" the local 
discrete solution (see (fT2]) ) [yoid, Ynew)'^ in the discrete problem defined 
by a higher order method. In such a way, one obtains an estimate r of 
the local truncation error. The same is done by considering an equivalent 
formulation of the same method (compare with (jl6p and (jl7p -()18p). thus 
obtaining a new approximation n = 7^"^ (8> Im'T- After that, the estimate 
e of the local error is obtained by performing one blended iteration, namely 
by formally solving the block diagonal linear system (see (1230 ) 

Ne = eT + {I -6)ti. 

In more detail, when using the GLM defined by the triple {k,r,£), the cor- 
responding method used for estimating r is the GBDF defined by the triple 
{k + l,r,i). Consequently, the cost for estimating e is the same as that for 
carrying out one blended iteration. This can be done for all methods listed 
in Tables [U and [21 with the only exception being the triple (3,2,2), for which 
the value of r should be increased. For this reason, such a method has not 
been considered for carrying out the numerical tests in Section [6l 

Finally, for the efficient variation of the stepsize, we have considered the 
Nordsiek implementation of the methods, firstly introduced in [2H] and later 
modified according to, e.g., \25\ I30j. For the sake of completeness, we also 
mention that the initial guess for the blended iteration (|22p - (|23p is obtained 
by extrapolation, through the interpolating polynomial on (see (fT2|l ) yoid- 

6 Numerical Tests 

In this section, we report a few numerical tests comparing a Matlab fixed- 
order implementation of the Blended GLMs presented here, with some of the 
most reliable codes currently available, on selected stiff test problems. Both 
the problems and the codes have been taken from the current release (re- 
lease 2.4) of Test Set for IVP Solvers [32]. In particular, we have considered 
the following solvers: 

• BiMD, which uses a blended iteration for solving the generated discrete 
problems; 

• GAMD, which uses a nonlinear splitting for solving the discrete prob- 
lems, generated by block BVMs in the Generalized Adams Methods 
(GAMs) family; 
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• DASSL, which implements standard BDF methods. 
The problems considered are: 

• Pollution (of dimension 20); 

• Elastic Beam (of dimension 80); 

• Emep (of dimension 66); 

• Ring Modulator (of dimension 15). 

The comparisons have been summarized in corresponding work-precision 
diagrams [32j, where the computational cost is plotted versus accuracy. A 
standardized cost has been computed as the number of floating-point oper- 
ations for the factorizations and the system solvings required by each code, 
while accuracy has been measured in terms of mixed-error significant correct 
digits |32j. For all codes, the used tolerances are essentially those specified 
in the Test Set. Figures [SHE] summarize the results obtained, where the la- 
bel "ORDER fe" is used for the fixed-order implementation of the Blended 
GBDF of order k (see Table [5]): as one can see, for each problem, the selected 
fixed-order implementations of the Blended GBDF presented here, appear 
to be competitive with the above mentioned solvers. 

7 Conclusions 

In this paper, a straightforward approach for deriving L-stable General Lin- 
ear Methods (GLMs) of arbitrarily high order has been presented. Such an 
approach relies on the class of Boundary Value Methods (BVMs) for ODEs. 
In the same framework corresponding starting procedures are easily derived, 
as well as appropriate error estimates. 

The generated discrete problems can be efficiently solved by means of 
the blended implementation of the methods, thus defining corresponding 
Blended GLMs, equivalent to the original methods, from the point of view of 
the stability and accuracy properties. The corresponding blended iterations 
are all L-convergent, thus appropriate for the underlying L-stable methods. 

A number of numerical tests, on problems taken from the Test Set for 
IVP Solvers, prove that the obtained methods are competitive with some of 
the best codes currently available. 

The availability of methods having arbitrarily high-order makes them 
good candidates for an efficient variable-order implementation. 
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Figure 2: Boundary loci of the GLMs obtained with the choice (jl4p . k 
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Figure 4: Boundary loci of the GLMs obtained with the choice (jlSp . k 
12,14,16. 
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Figure 5: Numerical results for the Pollution Problem. 
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Figure 6: Numerical results for the Elastic Beam Problem. 
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Figure 8: Numerical results for the Ring Modulator Problem. 
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